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ABSTRACT 

We study dark-matter halo density profiles in a high-resolution N-body simulation 
of a ACDM cosmology. Our statistical sample contains ^ 5000 haloes in the range 
IQ-'^-'^ — 10^'*/i~^Mq and the resolution allows a study of subhaloes inside host haloes. 
The profiles are parameterized by an NFW form with two parameters, an inner radius 
Ts and a virial radius i?vir, and we define the halo concentration Cvir = Rvit/^s- We 
find that, for a given halo mass, the redshift dependence of the median concentration 
is Cvir oc (1 + z)^^- This corresponds to rs{z) ~ constant, and is contrary to earher 
suspicions that Cvir does not vary much with redshift. The implications are that high- 
redshift galaxies are predicted to be more extended and dimmer than expected before. 
Second, we find that the scatter in halo profiles is large, with a la A(logCvii) = 
0.18 at a given mass, corresponding to a scatter in maximum rotation velocities of 
AKnax/Knax — 0.12. Wc discuss implications for modelling the TuUy-Fisher relation, 
which has a smaller reported intrinsic scatter. Third, subhaloes and haloes in dense 
environments tend to be more concentrated than isolated haloes, and show a larger 
scatter. These results suggest that Cvir is an essential parameter for the theory of 
galaxy modelling, and we briefly discuss implications for the universality of the TuUy- 
Fisher relation, the formation of low surface brightness galaxies, and the origin of the 
Hubble sequence. We present an improved analytic treatment of halo formation that 
fits the measured relations between halo parameters and their redshift dependence, 
and can thus serve semi-analytic studies of galaxy formation. 
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1 INTRODUCTION 

In the "standard" picture of galaxy formation, dark-matter 
(DM) haloes provide the framework for the formation of lu- 



minous galaxies (e.g., White fc Rees 1978 



Blumenthal et al 



1984; White fc Frenk 1991). The DM haloes are assumed 



to form hierarchically bottom-up via gravitational amplifi- 
cation of initial density fiuctuations. The haloes carry with 
them gas, which eventually cools and contracts to form lumi- 
nous disc galaxies at the halo centres. The halo profile has a 
direct dynamical role in determining the observable rotation 
curve of the disc. It also affects gas cooling and infall and 
therefore the structural properties of the resultant disc, such 
as size, luminosity and surface brightness. In order to model 
properly the dissipative stages of galaxy formation and ob- 
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tain meaningful predictions for observable quantities (such 
as the TuUy-Fisher relation), it is important to perform de- 
tailed dynamical studies of the evolution of halo structure, 
and to obtain statistical characteristics based on a fair sam- 
ple of the simulated halo population. 

Most naturally, the density profiles of haloes are ex- 
pected to be a two-parameter family. This is because, as- 
suming that the formation of haloes can be approximated 
by spherical collapse, each proto-halo perturbation can be 
characterized by two quantities, e.g., mass and radius (or 
density contrast) at some fiducial cosmological time. In the 
approximation of spherical collapse, these parameters spec- 
ify the full evolution of each halo, including the epoch at 
which it collapses and its virial radius. A successful two- 
parameter functional form for the halo profi l es ha s been 



proposed by Navarro, Frenk, & White (1995 
hereafter NFW95, NFW96, and NFW97): 



199e 



1997 
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pNFw('^) — 



{r/r,){l + r/r,y 



(1) 



where rs is a characteristic "inner" radius, and ps a corre- 
sponding inner density. As we show in § ^ one of the inner 
parameters can be replaced by a "virial" parameter, either 
the virial radius (iivir), mass (Mvir), or velocity (K-ir)- A 
very useful alternative is the concentration parameter Cvir, 
which relates the inner and virial parameters. NFW found 
that this functional form provides a good fit to haloes over a 
large range of masses, and for several different cosmological 
scenarios. It has been tested for the Einstein-deSitter model 
with a standard CDM power spectrum of initial fluctuations 
(SCDM), a flat cosmological model with flm ~ 0.3, JIa = 0.7 
and a corresponding CDM power spectrum (ACDM), and 
several mod els w ith power-law power spectra (confirm ed by 



Craig 1997, and Kravtsov, Klypin, & Khokhlov 1997, here- 



after KKK97). 

NFW then noticed that, for a given cosmology, their 
haloes show a strong correlation between the model's two 
parameters, e.g., an increase in pa for decreasing Mvir- A 
natural reason for the fact that low-mass haloes tend to 
show higher densities is that they typically collapsed earlier, 
when the universe was denser. To model this trend, NFW 
proposed a toy model (outlined in Appendix A) which as- 
sumes that Pa is a constant multiple k of the universal den- 
sity pu{zc) at a collapse redshift Zc, and that the collapsing 
mass at Zc is a constant fraction / of the total halo mass 
that has just collapsed and virialized at 2 = 0. The general 
trend of the relation between the two profile parameters at 
2 = is reproduced well for a proper choice of values for 
the constants k and /, with different values for the different 
cosmological models. 

Since the halo profiles are expected to be a two- 
parameter family, it is important to study the scatter about 
this mean relation between the two halo parameters. This 
scatter could provide the second parameter which is neces- 
sary in order to explain the observed variations in galaxy 
properties, such as bulge-to-disc ratio, size and surface 
brightness. It may be argued that the scatter in halo spin pa- 
rameter may also contribute to these variations, but it can- 
not account for all of them because, e.g., it does not properly 
correlate with the environment. The scatter in halo profiles 
should also have direct implications for understanding the 
surprisingly tight scatter observed in the TuUy-Fisher (TP) 
relation for disc galaxies. NFW found a small scatter among 
their simulated haloes, which could have provided a conve- 
nient explanation for the small TF scat ter, but other the- 
oret ical studies predict a larger scatter ( Eisenstein & Loeb 
1996), ; nd the generality of the NFW result is limited by the 
small number of haloes simulated per cosmology (~ 20) and 
by their selective choice of haloes near virial equilibrium. We 
therefore here study in detail the scatter in a large, "fair" 
sample of simulated haloes. A similar investigation has been 
performed by Jing (2000). 

The accumulating data of galaxies at high redshifts pro- 
vide a great incentive for studying the properties of the halo 
population as a function of redshift. NFW97 tried to extend 
their toy model in order to predict this redshift dependence 
by assuming that k and / are both constant in time. In 
order to actually study in detail the redshift dependence of 
halo profiles, we use our large statistical sample of simulated 
haloes. Our results below, which differ from the NFW97 toy- 



model predictions at 2: > 1, motivate modifications in the 
toy model in order to properly account for the simulated 
behavior. 

We present here an analysis of a statistical sample of 
halo profiles drawn from cosmological N-body simulations 
of ACDM with r2m = 0.3 and = 0.7. The unprecedented 
features of this analysis are: 

• The sample is large, containing about 5000 haloes in 
the mass range lO^-lO^^/i'^M© at 2 = 0. 

• The sampling is "fair", in the sense that haloes are 
found in any environment, field and clustered, and irrespec- 
tive of the dynamical stage of the halo after virialization. 

• The resolution is high, allowing a distinction between 
"distinct" haloes and "subhaloes", and a study of environ- 
mental trends. 

• The time evolution and scatter about the one- 
parameter family are studied in detail. 

In § ^ we discuss further the parametric functional form 
used for the halo profiles. In § |^ we present the revised toy 
model for predicting the mean relation between the halo 
profile parameters and its redshift dependence. In § ^ we 
describe our N-body simulations and our method of halo 
finding and classification, and discuss tests for the effects of 
mass resolution on fit parameters. In § ^ we present our re- 
sults for haloes at 2 = 0; we compare the mean result to our 
model prediction, and quantify the intrinsic scatter. In § ^ we 
discuss implications for observable rotation curves and the 
TF relation. In § |^ we investigate the redshift dependence 
of halo properties, and the toy-model fits. Finally, in § ^, we 
summarize our results and discuss further implications. 



2 PROFILE CHARACTERISTICS 

We choose to fit the density profiles of all haloes at all red- 
shifts with the NFW two-parameter functional form (Eq. 
This is a convenient way to parameterize the profiles, with- 
out implying that it necessarily provides the best possible 
fit. A similar analysis could be carried out using alterna- 
tive functional forms. In this section, we discuss the various 
parameters that are associated with the halo density pro- 
file, the relations between them, and how the values of these 
parameters influence observable quantities. 

The inner radius, rs, is where the effective logarithmic 
slope of the profile is —2, a characteristic radius which we 
term r_2. For much smaller radii, Pnfw oc r~^, and for much 
larger radii, Pnfw oc r~"^. The inner density parameter of 
the NFW profile is related to the NFW density at by 
Ps = 4pNFw(?'s), and equals the local density at about half 
rs: Ps = pNFw(r = 0.466 rs). 

The outer, virial radius -Rvir, of a halo of virial mass 
Mvir, is defined as the radius within which the mean density 
is Avir times the mean universal density pu at that redshift: Q 

Mvir = — AvirPuii?ir- (2) 



The associated virial velocity is defined byQ 



2 

vir 



t 

Mvir ^ lO"/i-iM0 (f7oAvir(2)/200)[ii:vir(l + z) /75 h-'^kpc]'^ . 
t Vvir ~ 75km/s(_Rvir/75 h-'^kpc){noA^i,{z)/200y/^ (1 + zf^^ . 
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Figure 1. Maximum velocity versus concentration. The maxi- 
mum rotation velocity for an NFW halo in units of the rotation 



GMvii/Rvir- The one-to-one relations between the three 
virial parameters are fully determined by the background 
cosmology (independent of the inner halo structure) , so only 
one of iViem at a time can .qerve in the pair of independent 
parameters characterizing the profile. The virial overdensity 
Avir is provid ed by the diss i pationless spherical top-h at col- 



lapse model (Peebles 1980 



Eke, Cole, & Frenk 1995); it is 



a function of the cosmological model, and it may vary with 
time. For the Einstein-deSitter cosmology, the familiar value 
is Avir — 178 at all times. For the family of flat cosmologies 
( ilm + ~ 1) , the val ue of Avir can be approximated by 
( [Bryan fc Norman 199^ ) Avir ^ (ISvr^ + 822; - 39x^)/n{z), 
where x = ^l{z) — 1, and i}{z) is the ratio of mean matter 
density to critical density at redshift z. For example, in the 

AriDM rngmnlngipal mnHpl that kptwk ag thp haaig fnr r>irr 



a.nalysi.'t in this paper (Q^ = n..S), the va.hie at z = is 
Avir ==; 337. 

An associated useful characteristic is the concentration 
parameter, Cvir, defined as the ratio between the virial and 
inner radii, 



Rvir/r 



(3) 



Note that our definition of Cvir differs slightly from that orig- 
inally used by NFW, Cnfw = R20o/rs, where i?200 is the 
radius corresponding to a density of 200 times the critical 
density, independent of the actual cosmological model. 

A third relation between the parameters of the NFW 
profile is 



Mvir = 47rpsr,^A(Cvir), A(Cvir) = ln(l + Cvir) - , .(4) 

1 + Cvir 

The three relations (Eqs. ^, | and ^ allow the usage of any 
pair out of the parameters defined so far (excluding degen- 
erate pairs consisting only of virial quantities) as the two 
independent parameters that fully characterize the profile. 

Since the more observable quantities have to do with 
rotation curves, it is worth translating the density profile 
into a circular velocity curve for the halo. 



= V' 

^ VI 



A(Cvir 



(5) 



where x = r/r^. The maximum velocity occurs at a radius 
Tmax — 2.16 and is given by 



1/2 

* m 



: 0.216 



A(Cvir 



(6) 



For typical Cvir values in the range 5 — 30, Knax varies in the 
range (1 — 1.6)Vvir. Figure 1 shows the ratio V^ax/V^ir as 
a function of Cvir. Note that for haloes of the same mass, a 
larger Cvir goes with a larger Vmax- Because of the relation- 
ship between rmax and rs, haloes with Cvir SlO have velocity 
curves which continue to rise gradually out to an apprecia- 
ble fraction of their virial radii, while those with Cvir ^10 rise 
more steeply, and possibly represent galaxies in which rmax 
is identifiable observationally (though the effects of baryonic 
contraction should be taken into account before the observed 
and simulated values of r-max can be compared). 

In order to gain a qualitative understanding of how pro- 
file characteristics may affect galactic disc formation, we 
may assume that an exponential disc forms by t he adiabatic 



al lis virial radluH tm a, fuiicLiuii of halo coiicuiiLrallou. 


et al. 


1986 




mers 


1997; 



Flores et al. 199 3 
1998^ 



Dalcanton, Spergel, & Sum- 



Mo et al. 



The final disc size, rj, can be 
derived from the halo parameters Cvir and i?vir, under the 
following further assumptions that ( a) the disc forms from 
cold gas of mass ~0.03A'/vir (see, e.g. Bomerville & Primack 



199!;) which follows the original density profile of the halo 
out to 7?vir, and (b) the specific angular momentum of the 
gas is equal to that of the halo, which has an original spin 
parameter of A = 0.035 — the most probable value of the 
well-known log normal distribution. Under these assump- 
tions, we find the foll owing fitting formula, is good to within 
1% for 1 < Cvir < 50 (iBullock et al. 2000b| ): 



5.7 h ^kpc 



i?vir 



100 fe-ikpc 



[1 + (Cvir/3.73)° 



(7) 



(A similar fitting formula, which allows more varied assump- 
tions abou t the halo and disc make-up, is provided by 



Mo 



et al. 1998 .) The general result is that rd, and thus the disc 
surface brightness, is a decreasing function of Cvir. 

We stress again that the specific choice of the NFW 
functional form does not limit the generality of our analysis 
in a severe way. This is largely due to the association of 
the specific with the more general r_2. When the NFW 
function is fitted to a generic halo whose profile even vaguely 
resembles a similar shape, the fitting procedure is likely to 
return an rs value that is close to the effective r_2 of that 
halo. The concentration parameter can then be interpreted 
as a general structure parameter not necessarily restricted 
to the specific NFW function. In particular, any spread in 
Cvir can be attributed to a real scatter in a "physical" inner 
radius such as r_2 rather than to inaccuracies in the assumed 
universal profile. 

The interpretation of rs as r_2 allows one to map 
the NFW parameters to appropriate parameters of other 
functional forms. For many purposes, such as determining 
Vinax or modelling gas cooling and galaxy formation, the 
NFW form is sufficiently accurate. However, a compari- 
son with alternative profiles with different core behaviors 
may be important when much smaller radii are concerned 
(r s0.02i?vir), where there are indications of deviations from 
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an extrapolation of the NFW shape ( Kravtsov et aL 1998 
" [Moore et al. 19 99 



hereafter KKBP98; 



Primack et al. 1999 



Klypin et al. 200C| , hereafter KKBPOO). 

A specific example of an alt ernative profile functional 
form is the Burkert profile (1995): 



PB(r) 



Pb 



[l + (r/rb)2][l + r/rb] 



(8) 



This profile resembles the NFW profile for r s,0.02_Rvir- The 
Burkert profile has a log slope of —2 at r_2 — 1.52 rb, so 
one can relate the scale radii of NFW and Burkert by rb — 
rs/1.52, and then relate the concentration parameters by 
Rvii/rb ~ 1.52 Cvir. 

Note also that the relationship between Cvir and rmax 
is robust, regardless of profile shape. For example, with the 
Burkert profile we have an implied velocity maximum at 
rmax — 3.25 rb. If one assumes instead the relation as gleaned 
from an NFW fit, with rmax — 2.16rs — 3.28rb, there is good 
agreement between the values of rmax. 



3 A REVISED TOY MODEL 

In any investigation based on computer simulations it is use- 
ful to have a simple toy model that helps interpret the nu- 
merical results and allows an easy application of the conclu- 
sions in subsequent analytic or semi-analytic investigations. 
As mentioned in the Introduction and outlined in Appendix 
A, NFW96 and NFW97 proposed such a model, with 2 free 
parameters, which successfully recovers the mean Cvir-Afvir 
relation at z = for the several different cosmological mod- 
els simulated by them. Our simulations (see §5 below) in- 
deed confirm the success of this model at 2: = 0. However, we 
find (§ ^ below) that the NFW97 model does not reproduce 
properly the redshift dependence of the halo profiles as seen 
in the simulation; it significantly over-predicts the concen- 
tration of haloes at early times, z^l. We therefore propose 
a revised toy model that is shown below to recover prop- 
erly the full behavior of the mean Cvir-Mvir relation and its 
redshift dependence. We present the model in this prepara- 
tory section, so that we can refer to it when we describe 
and interpret the results from the simulations in the follow- 
ing sections. A fortran code that implements this model for 
several cosmologies is available from the first author upon 
request, n 



3.1 The model 

We seek a model for the typical halo concentration, denoted 
in this section by Cvir(Mvir, a), for a given mass Mvir and 
epoch a — (1 + z)^^. Following the general spirit of the 
NFW97 model, we assign to each halo an epoch of col- 
lapse, ttc- Unlike their formulation, which utilizes the ex- 
tended Press-Schechter approximation, we define ac in a 
simpler way as the epoch at which the typical collapsing 
mass, M*(ac), equals a fixed fraction F of the halo mass at 
epoch a, 



M.(ac) = FMvir 



(9) 



The typical collapsing mass at an epoch a is defined by 
a[M*(o)] = 1.686/D(a), where a{M) is the a = 1 linear 
rms density fiuctuation on the comoving scale encompass- 
ing a mass M, D{a) is the linear growth rate, and 1.686 is 
the linear equivalent of the density at collapse according to 
the familiar spherical collapse model. The typical collapsing 
mass is therefore a known function of the linear power spec- 
trum of fluctuations and D{a) for the cosmology in hand. 
For some purposes we find it convenient to measure the halo 
mass in units of the typical halo mass at the same epoch. 



^ = A/vir(a)/M,(a). 



(10) 



The second relation of the model arises by associating 
density of the universe pu at with a characteristic density 
of the halo at a. NFW97 used the inner density parameter 
Ps, which is related to Mvir and rs via the specific shape of 
the NFW profile. Instead, we define a more general charac- 
teristic density ps by combining inner and virial quantities: 

Mvi. = ^r!p,. (11) 

The NFW profile implies pa — 3psyl(cvir). The association 
of the halo density pa with the universal density at collapse 
is made via a second free parameter, K: 

Ps = A'^Avir(a)p„(ac) = K'^ A^ir{a)pu{a) (^-^^ . (12) 

The parameter K represents contraction of the inner halo 
beyond that required for top-hat dissipationless halo viri- 
alization, and it is assumed to be the same for all haloes. 
Using Eqs. ^, |ll|, and we obtain a simple expression for 
Cvir in terms of ac as our second model relation: 

Cvir{p, a) = K — . 



(13) 



The model is fully determined by Eqs. ^ and |l3| given 
the values of the two parameters F and K. We find below 
(§ ^ and § ^) that by adjusting F and K we are able to 
reproduce the full behavior of Cvir(p, a) as measured in our 
simulations. The small differences in the definitions of F and 
K compared to the analogous parameters of NFW97, / and 
k, make a big difference in the success of the model. 

Note in Eq. ^ that, for any cosmology, ac is uniquely de- 
termined by Mvir, independent of a. This implies via Eq. ^ 
that, for a fixed halo mass, 



r(a) (X a. 



(14) 



Our model thus predicts that for haloes of the same mass the 
concentration is proportional to (1 -I- z)~^. This is different 
from the NFW prediction in which the concentration is a 
much weaker function of redshift. 

In order to gain a basic understanding of the impor- 
tant elements of this model, we discuss its predictions in the 
context of three cosmological models of increasing complex- 
ity: (1) a self-similar model of Einstein-deSitter cosmology 
and a power-law power spectrum of fluctuations, (2) stan- 
dard CDM, in which the universe is still Einstein-deSitter 
but the power spectrum deviates from a power law, and (3) 
the relevant cosmology of the current investigation, ACDM, 
with fim 7^ 1 and a non-power-law spectrum. 



http:/ /www. astronomy.ohio-state.edu/ james/CVIR/parts.html 
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3.2 Example 1: the self-similar case 

As an illustrative example, assume a fully self-similar case: 
Einstein-deSitter cosmology, f2m = 1, for which the growth 
rate is D{a) oc a, and a power-low power spectrum of 
fluctuations, P{k) oc fe", for which cr(M) oc M"" with 
Q = (n -f 3) /6. In this case, 

M.(a) oc a^/". (15) 

Together with Eqs. | and |l^, we have ttc/ao — {fiF)" . Then, 
using Eq. we obtain 

c,i,{li,a)=K{Ffi)-". (16) 

Note that in the self-similar case the two parameters F 
and K can be replaced by one parameter, KF~°'. Equiv- 
alently, we may vary only K and adopt the natural value 
F — 1, namely, apply the model to the collapse of the whole 
halo. This is a special feature of our revised model, not valid 
in the original NEW model. 

The slope of the function Cvir(M) at a is completely 
determined by the power index a (i.e., n): Cvir oc fJ.~°'. 
This simple mass dependence can be checked against the re- 
sults of the simulations of NFW97 in their Figure 6, which 
presents Cvir(M) at 2 = for Einstein-deSitter cosmology 
with four different power-law spectra: n = — 1.5, —1, —0.5, 
and (a — 1/4, 1/3, 5/11 and 1/2). In each case, our model 
predicts the simulated slope quite well, even slightly bet- 
ter than the NEW model, and, unlike the NEW model, it 
predicts an exact power law relation between Cvir and ^i, as 
would be expected for a scale-free cosmology. As the power- 
law M~" becomes steeper (n larger), the difference in col- 
lapse epochs for haloes of a given mass difference becomes 
larger, which is reflected in a steeper Cvir(M) relation. 

The collapse factor K may be determined from the sim- 
ulations at the present epoch by matching the value of Cvir at 
any desired /x. We find for a typical simulated halo {/j, = 1) 
at a = 1 that Cvir — 10, implying that the additional collapse 
factor for the whole halo {F = 1) must be A" ~ 10. 

Clearly, for the fully self-similar case one would expect 
any dimensionless properties of Af, haloes to be invariant 
in time. This is easily verified by setting = 1 in Eq. |l6[ 
Cvir(l,a) = K (for F = 1). In this case, the concentration 
is also fixed for any other fixed value of fi. Accordingly, the 
concentrations are different for haloes of the same mass that 
are addressed at different redshifts. For the Einstein-deSitter 
cosmology, Eq. ^ yields i?vir oc a, and then the general be- 
havior of Cvir oc a (Eq. |l^) implies (via Eq. |^ that the value 
of Ts is the same (in physical coordinates) at all redshifts. 
Again, this is different from the NFW97 model prediction. 

3.3 Example 2: SCDM 

As a second example, consider the standard CDM cosmology 
(SCDM: = 1, h ^ 0.5, erg = 0.7). Although the time 
evolution here is still self-similar, the power spectrum is not 
a power law: a{A'I) has a characteristic bend near ~ 1 
today [which corresponds to a mass M,{a = 1) ~ 3.5 x 
10" /i"^M0 ]. The local slope varies from a ~ for /i <^ 1 
to a = 2/3 for > 1. 

The model solution for Cvir(/i) does not have a closed 
form in this case, but it is easy to obtain a useful approx- 
imation. The slope of the Cvir(M) relation at a specific fi is 



determined by the effective slope of the power spectrum on a 
scale corresponding to the mass /iFM« (not fJ.Mt), since this 
is the mass scale used to characterize the halo collapse time 
m Eqs. |andjl^). We now obtain an approximate relation 
similar to Eq. but with the effective local a replacing the 
constant a: 

c.ir~J^(MF)-"(''^'. (17) 

Unlike the power-law example, the value chosen for the con- 
stant F does play a role in the slope of the Cvir{fi) relation at 
a given a. In order to determine the best values of F and K, 
we match the model predictions of Cvir(/i) to the results of 
the A''— body simulations of the SCDM model at the present 
epoch. Using F — 0.01 and K — 4.5 in Eqs. ^ and |l^, we 
are able to reproduce quite well the z — SCDM results 
of NFW97 (their Figure 6) over the range fi ~ 0.01 - 100. 
The relation about /i ~ 1 is well approximated using Eq. \vi{ 
where a(/iF = 0.01) ~ 0.14. 

But now, using no extra parameters, we are also able 
to reproduce the time dependence of this relation, which we 
test using the Cvir values as determined from our small box 
(7.5 h~^Mpc) SCDM simulation. In our simulation, we find 
consistent results with the NFW97 simulation at 2: = 0. In 
addition, for a fixed Mvir, we find that the model-predicted 
scaling of Cvir(a) oc a indeed describes very well the time 
evolution of the halo population. 

In our toy model, the values of F and K are assumed to 
be constants as a function of both a and jj.. Such a behavior 
is naturally expected in the fully self-similar case, in which 
no special time or scale is present in the problem. However, 
the success of this toy model in the SCDM case, which is 
not scale invariant, is somewhat surprising. The reason for 
this success is linked to the small value of F, which pushes 
the relevant mass scales of the problem to values much be- 
low that of the bend in the power spectrum (near ~ 1), 
where the spectrum is approaching a power-law behavior. 
For small values of fi (/isO.l), the slope of Cvir(A') is almost 
independent of the actual value of F as long as the latter is 
smaller than 0.05. The specific preferred value of -F = 0.01 
arises from the need to match the model Cvir(Ai) with the 
simulations in the range fi > 1. 

3.4 Example 3: ACDM 

Our third example concerns a currently popular ACDM cos- 
mological model (f2m = 0.3, flA ~ 0.7, h — 0.7, ug = 1.0, 
where A/* ~ 1.5 x 10"/i"^Mo at z ^ 0). In this case, self 
similarity is violated due to the non-power-law spectrum as 
before, and also by the time-dependent fluctuation growth 
rate associated with the low fim. Using similar reasoning to 
our discussion in the previous example, one may expect our 
toy model with constant F and K to break down. These 
worries are somewhat alleviated as long as F is small. As in 
the SCDM case, this pushes the relevant mass scales to the 
power-law regime away from the bend in the power spec- 
trum. In addition, a small value of F demands that the 
collapse epoch is early, when the mean density was near 
the critical value, nni(ac) — 1, and therefore the fluctuation 
growth rate was close to that in the self-similar cosmology, 
D{a) oc a. The deviation from the self-similar growth rate 
introduces in Eq. ^ a multiplicative correction factor, given 
by the growth of fluctuations between the epochs Uc and a 
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in the given cosmology compared to the Einstein-deSitter 
case. In the case of ACDM and Oc -C 1 this factor at a = 1 
is about 1.25. 

Although the small value of F alleviates most the ex- 
pected problems associated with the model, for open cos- 
mologies, there will be a break down at extremely high 
masses M 2> F~^Mt{a = 1). This is because the fluctuation 
growth rate saturates at late times, D(a ^ 1) ^ const., and 
the definition of collapse redshift (Eq. ^ loses meaning. For 
our best-fit value, F — 0.01, this limits the applicability of 
our model to haloes less massive than ~ 100M*(a = 1). For- 
tunately, haloes more massive than this are extremely rare 
(by definition) so the range of applicability covers almost all 
relevant cases. |^ 

As we show in § H and § |^ below, using F — 0.01 and 
K = 4.0 in the model equations ^ and ^ we are able to 
reproduce the full behavior of the median Cvir{fJ.,a) in our 
ACDM simulations. At a = 1, we have Cvii{fJ.) oc /i^™'^^', 
where a ~ 0.13 for jj. ^ 1. For haloes of fixed mass (pM, ~ 
const), we have Cvir(a) oc a as before. The same model pa- 
rameters also provide a reasonable fit to the z — ACDM 
resuhs of NFW97 over the range O.OIM. sMslOOA/.. 

We will show below that by setting the value of K to 2.6 
and 6.0 we are able to artificially parameterize the —la and 
+ la scatter respectively in the value of Cvir for the simulated 
population of haloes. A similar range of K values accounts 
for the scatter for all masses and at all cosmological epochs. 



4 SIMULATING HALOES 
4.1 The numerical simulations 

Only recently have large cosmological N-body simulations 
reached the stage where detailed structural properties of 
many dark-matter haloes can be resolved simultaneously. 
One of the most successful methods for high force resolu- 
tion and mass resolution is the Adaptive Refinement Tree 
(ART) code (KKK97). The method makes use of succes- 
sive refinements of the grid and time step in high density 
environments. The simulations based on the ART code pro- 
vide, for the first time, a compilation of a statistical sam- 
ple of well-resolved DM haloes, as well as substructure of 
haloes within haloes. In previous simulations, haloes were 
picked "by hand" using certain selection criteria from a low- 
resolution cosmological simulations, to be re-simulated with 
high resolution. This selection induces a certain bias into the 
sample. 

We have used the ART code to simulate the evolution of 
DM in a low-density flat ACDM model for which flm = 0.3, 
Qa ~ 0.7, h — 0.7, and as = 1.0 at 2 = 0. The simulation 
followed the trajectories of 256"^ cold dark matter particles 
within a cubic, periodic box of comoving size 60 h~^Mpc 
from redshift z = 40 to the present. We have used a 512'^ 
uniform grid, and up to six refinement levels in the regions 
of highest density, implying a dynamic range of 32, 768. The 

^ Choosing F = 0.001 and K = 3.0 extends the applicable range 
of our model to SlOOOA/»(a = 1), and results in only a slight 
over-prediction of Cvir values compared to the M ~ 10^"^ Mq 
haloes in our simulation, and a better fit to the ACDM haloes of 
NFW97. 



formal resolution of the simulation is thus /res = 1.8 /I'^kpc, 
and the mass per particle is nip = 1.1 x Mq . In the 

present paper, we analyze 12 saved time steps from z = 5 
to 0. We have also used two simulations in smaller boxes 
to check for resolution and cosmology dependence. One of 
these is a 30 h~^Mpc box simulation of the same ACDM 
cosmology, with 256^ particles, mp = 1.4 x Mq , and 

/res = 0.9 /i"^kpc. The other, in a 7.5 h ^Mpc box, is of 
the SCDM cosmology {flm = 1, ft = 0.5, and ag = 0.7 at 
2 = 0), and it consists of 128'^ particles, /res = 0.5 h'^kpc, 
and rup = 5.5 x lO^/i'^M© . Tests of the ART code for 
numerical effects on halo density profiles are discussed in 
KKBP98 and KKBPOO. 



4.2 Finding and fitting haloes and subhaloes 

In this investigation, we sample all types of DM haloes inde- 
pendent of their environment. In particular, we identify both 
the standard kind of "distinct" haloes, of the type identified 
using common friends-of-friends algorithms and considered 
in Press-Schechter approximations, and also "subhaloes", 
whose centres are located within the virial radius of a larger 
"host" halo. Our halo finding/classifying algorithm, which 
is based on the B ound Density Maxima technique (Klypin 
& Holtzman 1997), has been specifically designed to simul- 



taneously identify distinct haloes and subhaloes (Appendix 
B). 

We fit every identified DM halo using the NEW pro- 
file (Eq. ^). Before fitting, we check the halo radial density 
profile to see if it has a significant upturn, dp{r)/dr > 0, 
and if so, we declare this point to be the truncation radius 
Rt- Our measured Rt values are comparable or somewhat 
smaller than the expected tidal radii. For haloes with no 
significant upturn in density, we fit the NEW density pro- 
file using logarithmically spaced radial bins from 0.02i?vir 
to RviT, while for haloes with a truncation radius, we fit the 
profile from 0.02i?t to Rt, and extrapolate the NEW func- 
tion in order to assign fitted -Rvir and Mvir- 

The profile fitting is performed as follows. After identi- 
fying a centre for the halo, we count particles in each radial 
bin and assign corresponding Poisson errors based on the 
count in each bin. We then fit an NEW profile (by min- 
imization) to the counts in bins using the bin errors in the 
covariance matrix, and obtain best-fit values for the two free 
parameters iivir and Vs (or equivalently Mvir and Cvir, etc.) 
along with the corresponding errors in these parameters. We 
then remove unbound particles from each halo, as described 
in Appendix B, and iterate the process of determining Rt 
and fitting a profile until the halo contains only bound par- 
ticles. 

We present results for haloes with masses in the range 
1.5 X 10^^ — 10^'^h~^ Mq ; the smallest haloes thus contain 
^150 particles. A profile fit to a halo of only a few hun- 
dred particles may carry large errors. We therefore make a 
rigorous attempt to estimate the errors and take them into 
account in every step of the process. Poor fits due to small 
number statistics are marked by large errors that are incor- 
porated in the analysis and the results we present. 

The profile fit of haloes in crowded regions clearly in- 
volves ambiguities in the mass assignment to the subhaloes 
and the host. Our fitting procedure provides a well-defined 
recipe for mass assignment based on the value of Afvir 
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Figure 2. Convergence test for an Mvir = 2 X lO^^h-'^MQ halo 
simulated with 2, 000 and 120, 000 particles respectively. When 
the fit is restricted to 0.02iJvir — Rvir the best-fit Cvir values show 
no significant difference. 

even when the fit is actually performed inside an Rt that 
is smaller than the R^ir obtained by extrapolation.^] The 
concentration parameter is defined in the same way for all 
haloes, Cvir ~ Rvir/rs- Because in most cases of subhaloes the 
extrapolation procedure adds much less mass than the mass 
that actually lies between Rt and i?vir, the double counting 
is not severe; most of the mass associated with the upturn in 
the profile is assigned to a different subhalo or to the host. 
On the other hand, a small subhalo does not cause a sig- 
nificant upturn in the profile of its host halo, and its mass 
is therefore also included in the mass assigned to the host. 
This partial double counting introduces some uncertainty to 
any recipe for assigning a luminous galaxy mass to a halo of 
a given mass. 

The outcome of the halo finder/classifier is a statistical 
halo catalog that includes all the bound virialized systems 
in the simulation above the minimum mass threshold. We 
include distinct haloes and subhaloes, but not subhaloes of 
a second generation, i.e., those whose hosts are themselves 
subhaloes of a larger host. The output for each halo includes 
the list of its bound particles, the location of its centre, 
the NFW profile parameters (e.g., Cvir and Mvir), the corre- 
sponding errors (ctc and ctm), and the truncation radius, if 
relevant. 

4.3 Tests of resolution 

In this section we discuss the efi'ect of varying mass resolu- 
tion on the measured fit parameters of haloes. Resolving the 
detailed shape of halo density profiles at radii r s0.01i?vir 
requires a sufficient number of particles within the halo 

II For 5% of the subhaloes we actually find Rt < rg. In these cases 
the errors in the extrapolated values of iJvir and Mvir become 
especially large, but they are properly taken into account in the 
analysis. 
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Figure 3. Convergence test for Cvir evolution and scatter. Shown 
is a comparison of Mvir = 3 — 10 X 10^^ /i~^Mq haloes simulated 
using our main simulation (thick lines) and a second simulation 
with 8 times the mass resolution (thin lines). The solid lines and 
errors reflect the median and Poisson uncertainty respectively. 
The dashed lines reflect the estimated intrinsic scatter. There 
is no evidence for significant deviations in either the measured 
median or scatter as the mass resolution is increased. 



(KKBPOO; Moore et al. 1999), and one may worry that 
our mass resolution will inhibit our ability to determine 
density profiles for all but the most massive halos in our 
sample. However, a study of the inner slope of density pro- 
files at small radii is not the aim of this investigation. We 
aim to characterize the general shape of halo profiles from 
0.02i?vir — Rvir by measuriug Cvir values of halos, and our 
tests explore how mass resolution affects results when the 
fits are constrained to this outer radial range. 

The first test of mass resolution effects was performed 
using a multiple mass resolution version of the ART code 
(KKBPOO) in which a halo of interest is simulated with fine 
mass resolution, while the surrounding regions are simulated 
with particles of masses that are increasing as a function of 
distance from the halo. The simulation of a particular halo 
was performed with varying maximum mass resolution. In 
Figure |^ we compare density profiles of the same halo of 
mass 1.5 x 10^^h~^ Mq at z = in the highest and the 
lowest resolution runs that simulated a box of 30ft~^Mpc on 
a side. The cosmology was a low-density flat ACDM model 
for which Sim = 0.3, Qa = 0.7, h = 0.7, and ag — 0.9. 

The high mass resolution simulation was run with a 
peak force resolution of 0.2 h~^kpc and 500,000 timesteps 
at the deepest refinement level. The smallest mass parti- 
cle in this run has a mass of 1.7 x 10'^ Mq , implying 
an effective particle number of 120, 000 within the virial ra- 
dius of the halo. The setup was such that no particles other 
than those of the smallest mass ended up within the cen- 
tral ~ 100 ft~^kpc. The lowest resolution run was similar 
except that the mass resolution was such that the final halo 
had 60 times fewer (2000) particles within its virial radius. 
We fitted the resultant halo profiles in the two runs using 
the same procedure outlined in the previous section and ob- 
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tained Cvir = 17.45 ± 0.05 and Cvir = 17.70 ± 0.70 for the 
highest and the lowest mass resolution runs, respectively. 
The fits are shown in the two panels of Figure GI The ob- 
tained values of Cvir are consistent within errors, despite the 
vastly different mass resolutions. This gives us confidence in 
our fitting procedure and resulting values of the concentra- 
tion parameter. 

In order to extend our tests to smaller particle numbers 
and to test the dependence of the measured scatter and red- 
shift behavior of fitted Cvir values (see §5.3 and 7), we have 
compared the results of the ACDM simulation used for our 
main analysis with a simulation with the same cosmologi- 
cal parameters but with half the box size and 8 times the 
mass resolution. This simulation was stopped at 2 = 1.7, 
so we pursued this test only at that epoch and earlier. Be- 
cause the simulation of higher resolution contains only few 
haloes near the high mass end, we limit the comparison to 
the mass range (3-10)x10"/i"^Mq . The implied particle 
numbers for the low and high resolution simulations are 
300 - 1,000 and 2,400 - 8,000 respectively. The evolution 
of measured Cvir values as a function of z for the two simu- 
lations is shown in Figure ^. The solid lines are the median 
relations, the error bars reflect the Poisson uncertainty from 
the finite number of haloes, and the dashed lines are the 
estimated intrinsic scatter (see §5.3 for a discussion of our 
determination of intrinsic scatter). In each case, the thick 
lines correspond to the lower mass resolution simulation and 
the thin hues correspond to the higher resolution case. Both 
the median and the intrinsic spread in Cvir agree to within 
~ 5%. This level of agreement is consistent with the Pois- 
son error of the high-resolution simulation, which is of order 
10% in the mass range we consider. 

We conclude that our mass resolution is perfectly ade- 
quate for the purposes of this paper. 



5 HALO PROFILES TODAY 

We start by studying the halo profiles at the current epoch 
in the simulation. First, we study the median Cvir (which 
is also close to its mean) as a function of mass. Then, the 
dependence on the environment and on being a subhalo is 
presented. Finally, we discuss the scatter about the median 
Cvir(M), which leads to a discussion of possible implications 
on the observed rotation curves and the Tully-Fisher relation 
in the following section. 

5.1 Median relations for distinct haloes 

Figure ^ shows Cvir(Afvir) at z = for distinct haloes. The 
Poisson errors reflect the number of haloes within each mass 
bin. In order to account for the fit errors, we generated 500 
Monte Carlo realizations in which the measured Cvir and 
Mvir of each halo were perturbed using random Gaussian 
deviates with standard deviations equal to Uc and itm re- 
spectively. Median values of Cvir were then determined using 
the Monte Carlo ensemble. The lowest-mass and highest- 
mass bin have ~ 2000 and 20 haloes respectively (we avoid 
using bins with fewer than 10 haloes.), and the Poisson er- 
rors grow with mass accordingly. 

The median Cvir decreases with growing mass, in qual- 
itative agreement with the toy models, and therefore con- 
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Figure 4. Concentration versus mass for distinct haloes at 2 = 0. 
The thick solid curve is the median at a given Mvir. The error 
bars represent Poisson errors of the mean due to the sampling of a 
finite number of haloes per mass bin. The outer dot-dashed curves 
encompass 68% of the Cvir values as measured in the simulations. 
The inner dashed curves represent only the true, intrinsic scatter 
in Cvir, after eliminating both the Poisson scatter and the scatter 
due to errors in the individual profile fits due, for example, to the 
finite number of particles per halo. The central and outer thin 
solid curves are the predictions for the median and 68% values by 
the toy model outlined in the text, for F = 0.01 and three different 
values of K. The thin dot-dashed line shows the prediction of the 
toy model of NFW97 for / = 0.01 and fc = 3.4 X 10^. 



sistent with the assertion that small mass haloes are more 
concentrated because they typically collapse earlier than 
haloes of larger masses, rj The NFW97 model outlined in 
Appendix A has been slightly adjusted to yield Cvir rather 
than Cnfw. Its predicted slope is in reasonable agreement 
with that derived from the simulations, but there is some 
indication that the slope is too shallow. Using our revised 
toy model outlined in § ||, with F = 0.01 and K = 4.0, 
we reproduce the median relation even better, as shown 
by the central thin solid line. Near /i ~ 1 (Mvir ~ M, ~ 
1.5 X 10^^h~^ Mq at z = 0), the model prediction is 

Cvir(/i, 2 = 0) ~ 1.25K (^iF)-°(^-^' ~ 9Ai-°-". (18) 

Indeed, the slope of the Cvir(/i) curve is closely related to the 
varying slope of the mass power spectrum, which influences 
the relative difference in collapse epochs for typical objects 
on different mass scales. The factor of 1.25, as explained in 
§ ^ is a measure of the deviation from the Einstein-deSitter 
self-similar linear ffuctuation growth rate D{a) oc a between 
some high redshift and z = (where the corresponding col- 
lapse epoch, for a given mass, is earlier in the ACDM case). 



** The possibility that this trend is accentuated by the tendency 
for substructure to be better resolved in large-mass haloes has 
been investigated. We conclude that this is not the case, since 
we observe that the number of subhaloes within fixed mass hosts 
does not correlate with the measured concentrations. 
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Figure 5. Concentration versus mass for sublialoes at 2: = 0. Tlie 
curves and errors are the same as in Figure jij 



5.2 Subhaloes and environmental dependence 

If the median Cvir(Afvir) indeed reflects different formation 
epochs, one might expect the Cvir of haloes of a given mass 
to vary with the density of the environment, since haloes 
in dense regions typically collapse earlier. In particular, the 
concentration should tend to be larger for subhaloes com- 
pared to distinct haloes. Another effect that may lead to 
higher Cvir values in subhaloes is the exp ected steepening of 
their outer profile due to tidal stripping ( Klypin et al. 19991 



Ghigna et al. 199^ ; Okamoto fc Habe 199^ ; Avila-Reese et al 



199!:). Since stripping is likely to be more effective for small 
mass haloes, this process may lead to a stronger mass de- 
pendence in subhaloes. A third effect is that haloes that are 
embedded in a high-density environment are likely to expe- 
rience extreme collapse histories and frequent merger events 
which may affect their final concentrations. 

Figure ^ shows the relation Cvir(Afvir) at z = for sub- 
haloes. We see that subhaloes on galactic scales (Afvir ~ 
are indeed more concentrated than distinct 
haloes of the same mass. The dependence on mass seems 
to be stronger for subhaloes than for distinct haloes, with 
Cvir oc /i"'^ '^ (compared to fj,~"'^^), though the large errors 
in the case of subhaloes make this trend only marginal. 

We address directly the dependence of concentration 
on background density in Figure ^, which shows Cvir as a 
function of local density for all haloes (both distinct and 
subhaloes) in the mass range 

The local background density is defined for the dark mat- 
ter within spheres of radius l/i~^Mpc in units of the aver- 
age density of the universe in the simulation, pu = 8.3 x 
IO^^/i^Mq /Mpc^ We see that haloes in more dense envi- 
ronments indeed tend to be more concentrated than their 
more isolated counterparts. Note that this trend is in fact 
stronger than the dependence of Cvir on mass. 

We find a similar trend when the local density is de- 
termined within spheres of radius 1.5/i~^Mpc, but the trend 
becomes weaker when spheres of radius 0.5/i~^Mpc are used. 
Similarly, we find that the trend holds for haloes of mass 
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Figure 6. Concentrations versus environment. The concentration 
at 2 = of all haloes in the mass range 0.5 — 1.0 X lO^^/i~^M0 
as a function of local density in units of the average density of 
the universe. The local density was determined within spheres of 
radius l/i~^Mpc. The solid line represents the median Cvir value, 
the error bars are Poisson based on the number of haloes, and the 
dashed line indicates our best estimate of the intrinsic scatter. 



s5 X 10^^/i~^A'/q . For larger masses, the trend seems to 
become less pronounced, but this is quite inconclusive be- 
cause we have only a few massive haloes (sl5) with local 
densities s^lOOpu. A similar trend has been seen for distinct 
haloes alone, but it is difficult to make a definitive assess- 
ment in this case because there are only a few distinct haloes 
with local densities alOOpu. 

The results at low and high densities are comparable to 
those for distinct haloes and subhaloes respectively, consis- 
tent with a significant correlation between being a subhalo 
(or a distinct halo) and residing in a low-density (or a high- 
density) environment. 

Our toy model is not sophisticated enough to account 
for the dependence of concentration on environment, and for 
the exact relation for subhaloes. However, in the framework 
of our toy model (§ |^, we can artificially parameterize these 
trends by varying the collapse parameter A" as a function of 
local density. We are left with the qualitative speculations 
mentioned above for the interpretation of the trends with 
environment seen in the simulations. 

We discuss possible implications of the environment de- 
pendence in § H. 



5.3 Scatter in the concentration parameter 

One of the most interesting of our results at z = is 
the spread in Cvir values for fixed Afvir. A significant scat- 
ter about the median relations may have intriguing obser- 
vational implications. Before we report on our results, we 
briefiy describe our methods of ascertaining the intrinsic 
scatter. 

There are two sources of scatter on top of the intrinsic 
spread in halo concentrations. First, Poisson noise due to the 
sampling of a finite number of haloes in each mass bin adds a 
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significant spread, especially in the case of large-mass haloes. 
Small-mass haloes, on the other hand, are plentiful, but the 
relatively small number of particles in each halo introduces 
a significant error in the measured value of the halo profile 
parameters Mvir and Cvir- This is the second source of addi- 
tional scatter. The Poisson error due to the finite number of 
haloes is relatively straightforward to correct, but correcting 
the error in the profile parameters requires a more involved 
procedure. 

We account for this error in the profile parameters us- 
ing the errors obtained in the profile fits. Within each mass 
bin, we have performed ~ 500 Monte Carlo realizations in 
an attempt to undo the effect of the measurement errors 
as follows. Every measured Cvir value has been perturbed 
by a one-sided random Gaussian deviate, positive or nega- 
tive depending on whether the measured Cvir is smaller or 
larger than the median respectively. The standard deviation 
of each Gaussian deviate was set to be the error in the value 
of Cvir as estimated in the profile fit of that specific halo. 
The smaller scatter obtained in this set of Monte Carlo re- 
alizations provides an estimate for the spread excluding the 
fit errors. We then subtract in quadrature the Poisson error 
due to the finite number of haloes to obtain our estimate for 
the intrinsic scatter in Cvir. 

We have checked our technique of measuring the intrin- 
sic spread using an artificial ensemble of 1000 (spherical) 
haloes with a variety of numbers of particles and a known 
intrinsic distribution of Cvir. The technique reproduced the 
median concentration and true spread to within 5% when 
the particle number was varied from 100 to 10^, the range 
of interest for our simulated haloes. 

As discussed in §4.3, we have also checked our proce- 
dure for measuring the intrinsic scatter using a simulation 
of higher resolution in a smaller box of side 30 fe~^Mpc, 
in which there are on average 8 times as many particles 
in a halo of a given mass. Although the measured scatter 
was larger in the lower resolution simulation, the estimated 
intrinsic intrinsic spreads in the two simulations agree to 
within ~ 5% (see Figure |^ . This gives us confidence in our 
technique. 

Note that what we treat as measurement error in the 
profile fit actually includes scatter due to real deviations 
of the halo structure from a purely spherical NEW profile, 
which we should probably regard as part of the intrinsic 
scatter. This means that our estimated intrinsic scatter is a 
conservative underestimate. 

Now that the method has been discussed, we turn the 
attention back to the relation between Cvir and Mvir for dis- 
tinct haloes, Eigure ^. The measured 68% scatter is shown, 
as well as the "pushed in" corrected scatter which marks 
our (under-) estimated intrinsic scatter. As can be seen by 
noting the Poisson error bars, the correction at the small- 
mass end is almost entirely due to the measurement errors 
of the profile parameters, which are dominated by the small 
number of particles per halo. 

We see that the intrinsic spread is large; it is compa- 
rable to the systematic change in the median value of Cvir 
across the entire mass range studied. In addition, the spread 
is roughly constant as a function of mass, with a la devia- 
tion of A(logCvir) ~ 0.18. We discuss possible observational 
implications of this scatter in the next section. 

The spread in Cvir values as a function of Afvir for sub- 
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Figure 7. The probability distributions of distinct haloes (solid 
line) and subhaloes (dashed line) at z = within the mass range 
(0.5 - 1.0) X lO^^h-^MQ . The simulated distributions (thick 
lines) include, the ~ 2, 000 distinct haloes and ~ 200 subhaloes 
within this mass range. Log-normal distributions with the same 
median and standard deviation as the measured distributions are 
shown (thin lines). Subhaloes are, on average, more concentrated 
than distinct haloes and they show a larger spread. 

haloes is shown in Eigure ^. Note that the scatter is larger for 
the subhalo population than for their distinct counterparts 
of the same mass, with a la variation of A(logCvir) ~ 0.24. 
This is clearly seen in Eigure |^ where the probability distri- 
butions of concentrations for distinct haloes and subhaloes 
are compared (for Mvir = 0.5 - 1.0 x 10^^/i"^Mq ). It is 
possible that the larger scatter evaluated for subhaloes is a 
result of their more complicated formation histories, includ- 
ing for example more interactions and stripping. We point 
out that we have found no significant trend with the number 
of co-subhaloes within the same virialized host. Such a trend 
might have been expected if interaction among co-subhaloes 
plays an important role in determining the profile shape. 

The 68% intrinsic spread in Cvir as a function of the lo- 
cal density (for 0.5 - 1.0 x 10^^/i"^Mq haloes) IS given m 
Eigure ^ We can use the obtained distribution of halo con- 
centrations as a function of local density to probe questions 
associated with the origin of LSB galaxies and the observed 
morphology density relation (§ ^). 

In Eigure ^ we show the distribution of concentration 
values for distinct haloes and subhaloes, along with log- 
normal functions with the same median and standard de- 
viation. The log-normal forms describe the observed distri- 
butions reasonably well. Such a result has also been reported 
by Jing (2000). Our measured scatter for distinct haloes is 
similar to that reported by Jing for equilibrium haloes. 

In the context of the toy model presented in § ^, one can 
parameterize the spread in Cvir as spread in collapse epochs 
and/or collapse histories, via the parameters and K re- 
spectively. Using Eq. we find that the evaluated spread 
in Cvir can be matched by a spread of A[log K {ao / a^)] — 
A(logCvir) ~ 0.18 in the toy model. If we absorb all the 
scatter in the collapse parameter K, we find that the model 
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Figure 8. The spread in NFW rotation curves corresponding 
to the spread in concentration parameters for distinct haloes of 
3 X 10^^ Mq at 2: = 0. Shown are the median (solid), itlcr 
(long dashed), and ±2(t (dot-dashed) curves. The corresponding 
median rotation curve for subhaloes is comparable to the upper 
Icr curve of distinct haloes. 



matches the 50 ± 34% (encompassing 68%) percentiles of 
the Cvir distribution with K = 6.0 and 2.6 respectively (for 
F = 0.01). Note that the scatter in K is not symmetric 
about the median (of K — 4.0); it rather reflects the log- 
normal nature of the Cvir distribution. The model predic- 
tions for the above values of K are shown as thin solid lines 
in Figure ^ they match the simulated scatter fairly well. 
We show in § 1^ that this parameterization also reproduces 
the simulated scatter as a function of z. This is just a use- 
ful parameterization of the scatter using the toy model. A 
more detailed modelling of the spread with deeper physical 
insight is beyond the scope of the present paper. 



6 ROTATION CURVES AND TULLY-FISHER 

The simulated distributions of Cvir values as a function of 
mass and environment have several observational implica- 
tions. Here, we discuss only preliminary predictions involv- 
ing rotation curves and the TuUy-Fisher relation based on 
very crude assumptions about associating disc properties to 
those of the simulated dark-matter haloes. A more detailed 
study requires realistic modelling of physical processes in- 
volving gas and stars. 

In order to illustrate what the spread in Cvir values may 
imply observationally. Figure ^ shows example NFW rota- 
tion curves for 3 x 10^^ Mq (distinct) haloes using the 
median, ±la, and ±2(t values of Cvir. These are raw rotation 
curves of the dark-matter haloes before they are affected by 
the infall of baryons, but they may still serve as a crude ap- 
proximation for the final rotation curves. One can see that 
the rotation curves span a significant range of shapes and 
the corresponding spread in Vmax values is substantial. The 
median rotation curve for 3 x 10^^h~^MQ subhaloes (not 
shown) is similar to the upper la curve shown in Figure H 



A clue for the expected TF relation of discs may be pro- 
vided by the measured relation between the halo parameters 
Mvir and Vmax. The latter is derived from Cvir using Eq. ^ 
The A/vir-Knax relation is shown in Figure ^ separately for 
distinct haloes and subhaloes. The median distinct-halo re- 
lation is well approximated by the linear relation 



log[Mvir/(/i"'M0 )] = alog[V;„ax/(km/s)] +/3 



(19) 



with Of ~ 3.4 ± 0.05 and /3 ~ 4.3 ± 0.2, where the Poisson 
errors in each mass bin have been propagated to obtain the 
quoted error on each fit value. Note that the slope is steeper 
than that expected from the standard scaling of the virial 
parameters: Mvir oc V^i,.. This is a direct result of the cor- 
relation between mass and concentration (Mvir oc c~^'^^). 
We may in fact derive the expected a using the effective 
power law from Eq.^: Knax/K-ir oc Cvir^; which implies 

Mvir (X K^,„(Kir/Knax)' « V^^^C^^''^ OC F^,!. We poiut 

out that the linear relation provides a good fit, showing no 
obvious need for non-linear corrections in the TF relation. 

The lower panel in Figure ^ shows Mvir vs. Vmax for 
subhaloes. This relation is also well fit by the linear relation, 
Eq. |l|, but now with q ~ 3.9 ± 0.25 and /3 ~ 2.6 ± 0.75. The 
subhalo relation has a steeper slope compared to distinct 
haloes, and Mvir — W^^ A4q subhaloes typically have a 
~ 12% higher Knax. This difference between the slope and 
zero-point of distinct haloes and subhaloes may have impli- 
cations for the use of cluster or group galaxies to calibrate 
the TuUy-Fisher relation in the field. 

We point out, however, that if for subhaloes we replace 
Mvir by the mass Mt inside the truncation radius Rt, the 
logarithmic slope becomes at = 3.6 ± 0.2, consistent with 
the slope obtained by Avila- Reese et al. (1999) for haloes 
within clusters using a similar mass assignment procedure. 
(The reason for the slope change is that the ratio of Mt /Mvir 
is roughly 1 for low-mass, high-Cvir haloes, and becomes less 
than 1 for high-mass, low-Cvir haloes.) This slope is similar 
to the slope we find for distinct haloes. It is not obvious a 
priori which of the halo masses is more relevant to the mass 
of the cooled gas that ends up as the luminous disc, and 
thus to the observed TuUy-Fisher relation. Therefore, the 
worry about the universality of the slope of Tully-Fisher is 
not conclusive. However, the zero-point difference between 
the two types of haloes exists regardless of the mass choice, 
and is a robust result. 

The scatter in the TF relation is an issue of great inter- 
est. We find for distinct haloes a la scatter at fixed Vmax of 
A (log A/vir) — 0.17, while the corresponding scatter for fixed 
A/vir is AVinax/Knax — 0.12. This scattcr is in rough agree- 
ment with the spread predicted by Eisenstein & Loeb (1996) 
for a similar cosmology using Monte Carlo realizations of 
halo formation histories based on the Press-Schechter ap- 
proximation. The subhalo relation shows an even larger scat- 
ter, with AVinax/Knax — 0.16 at a fixed Mvir. 

Observational estimates for the in trinsic scatter in (I - 
band) TF r ange from a(V) /V ~ 0.09 ( [Willick et al. 199e| ) 
to ^ 0.03 (Bernstein 1994). At best, the observed scatter 



leaves no room for any intrinsic variation in the mass-to- 
light ratio of galaxies, and may imply that gas contraction 
and other hydrodynamical processes must somehow act to 
decrease the scatter. A simple idea that may resolve this 
discrepancy is discussed in Appendix C and we briefly out- 
line the argument here. For a flxed mass and spin, a more 
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Figure 9. Virial mass versus Vmax at z = for distinct lialoes 
(a) and subhaloes (b). The outer dot-dashed and dark dashed 
lines indicate the measured and corrected intrinsic 68% scatter 
respectively. 
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Figure 10. Median Cvir values as a function of Mvir for distinct 
haloes at various redshifts. The error bars are the Poisson errors 
due to the finite number of haloes in each mass bin. The thin 
solid lines show our toy model predictions. 



concentrated halo (higher Knax) will induce more gas con- 
traction, and therefore produce a smaller, brighter disc. Such 
a correlation between the mass-to-light ratio of galaxies and 
the deviation of Knax from the median Vmax at a given Mvir 
could reduce the scatter to that required to match observa- 
tions. Detailed modelling, including the back-reaction of the 
halo during disc formation, is needed to test this hypothesis 
in detail. 



7 REDSHIFT DEPENDENCE 

As data accumulate at high redshift, it becomes increasingly 
important to study the predicted evolution of the population 
of halo profiles as a function of redshift. 

Figure ^ shows the median Cvir as a function of Mvir 
for distinct haloes at several different redshifts. We see that 
for a fixed mass, the typical Cvir value changes quite dra- 
matically, while the shape of the mass dependence remains 
roughly constant. The thin solid lines show our toy model 
predictions. This two-parameter model, which has been nor- 
malized to match the slope and normalization of the relation 
at 2; = 0, does remarkably well at all redshifts. As predicted 
by the toy model in § bl the concentration of haloes of a 
fixed mass scales (1 -t z) A similar behavior has 

been confirmed using the SCDM simulation in a smaller box 
(7.5/i~^Mpc, described in § |^. The redshift dependence of 
the subhalo concentrations seems similar, but we don't have 
sufficient statistics for conclusive results involving subhaloes 
at high redshifts. 

As mentioned in the Introduction, the dramatic evolu- 
tion in the concentration of haloes of a fixed mass is differ- 
ent from the prediction of the NFW97 analytic toy model 
(see Appendix A). This is illustrated in Figure |ll|, which 
shows the median Cvir of the distinct halo population of 
Mvir = (0.5 - 1.0) X IQ^'^h-^MQ as a function of redshift. 
The NFW prediction (for 0.8 x lO^^/i'^M© haloes) overes- 



timates Cvir by ~ 50% aX z = 1, and the disagreement grows 
with redshift. p| Since NFW's approach was to resimulate a 
small number of halos identified at z = in a larger simu- 
lation, they could not check the redshift dependence of the 
halo concentration. Their resimulations typically had not yet 
collapsed to a single halo at high redshift. Our revised toy 
model reproduces the simulated redshift trend very well. The 
scatter about the relation is remarkably constant as a func- 
tion of z: A(logCvir) ~ 0.18. Also shown is how the spread 
can be parameterized by varying K in our toy model, as 
discussed in § ^ 

The redshift dependence of the inner radius, r^, can be 
deduced from that of Cvir by recalling that the virial radius 
of fixed-mass haloes also varies like J?vir oc aV^^/(1 + z). 
This implies that, on average, the inner radius of haloes of a 
given mass remains roughly constant as a function of redshift 
(aside from the z dependence of Avir). We see this explicitly 
in Figure which shows the evolution of the median and 
68% scatter of as a function of z for distinct haloes in 
the mass range 0.5 - 1.0 x lO^^/i"^M0 . The fact that the 
median value declines slightly near 2 = is due to the 
z dependence of Avir in the ACDM model simulated. The 
robustness of the characteristic length scale, r^, may provide 
an interesting clue for the understanding of the build-up of 
DM halo structure. 

The strong decline in the concentration of haloes of a 
fixed mass as a function of redshift should have an inter- 
esting impact on galaxy-formation modelling at high red- 
shift — e.g., ai med at understanding the nature of Lyman 
Break Galaxies (steidel et al. 1996; Lowenthal et al. 1997) 



tt As discussed in Appendix A, if the NFW model is modified 
by setting the free parameter / to the unphysically small value 
/ = 10~^^, then the redshift behavior can be reconciled with 
what is observed. However, for this small value of /, the NFW 
model predicts an excessively flat Cvir vs. Mvir trend. 
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Figure 11. Concentration as a function of redshift for distinct 
haloes of a fixed mass, Mvir = 0.5 - 1.0 X lO^^h'^M© . The 
median (heavy solid line) and intrinsic 68% spread (dashed line) 
are shown. The behavior predicted by the NFW97 toy model is 
marked. Our revised toy model for the median and spread for 
8 X 10^^ Mq haloes (thin solid lines) reproduces the observed 
behavior rather well. 



8 CONCLUSIONS AND DISCUSSION 

The main direct conclusions of this paper, based on ana- 
lyzing a statistical sample of dark-matter haloes in a high- 
resolution simulation of the ACDM cosmology, are as fol- 
lows: 

• The redshift dependence of the halo profile parameters 
has been measured in the simulations, and reproduced by an 
improved toy model. For example, Cvir oc {l + z)~^ for haloes 
of the same mass, predicting that at high redshift they are 
less concentrated and with larger inner radii than previously 
expected. The corresponding prediction for rotation curves 
is lower values of Vmax/Vvir at high z. 

• The correlation between any two halo profile parame- 
ters has a significant scatter. For example, in the Cvir-A/vir 
relation, the spread in Cvir is comparable to the systematic 
change in Cvir across three orders of magnitude in Mvir. The 
la spread for fixed Mvir is A(logCvir) — 0.18, corresponding 

to AVmax/Vinax — 0.12 at a givCU Mvir- 

• There are indications for environmental trends in halo 
profiles. Haloes in dense environments, or subhaloes, are 
more concentrated than their isolated counterparts of the 
same virial mass, and they exhibit a larger scatter in Cvir. 

The main implications of the above findings can be sum- 
marized as follows: 



and the evolution of the TuUy-Fisher relation ( Vogt et al 
1997|)7^1though, in general, haloes, and therefore galaxies, 
are expected to be smaller at high redshift (refiecting the 
higher universal density) and to have higher circular veloc- 
ities (Vvir oc R'^y^), the observed Cvir behavior will tend to 
counteract this tendency. 

Insight into the expected TF evolution of discs may be 
gained from Figure |l^, which shows the Mvir versus V^max re- 
lation for (distinct) haloes at several redshift steps. The evo- 
lution in the zero-point is indeed less dramatic than would 
be expected from the scaling of Vvir. In fact (not shown) 
there is almost no evolution in the zero-point between z = 
and 0.5. The slope of the relation is roughly constant as a 
function of redshift (a — 3.4±0.1) and the scatter is roughly 

constant; AVmax/Vmax — 0.12. 

Furthermore, because disc size is expected to be a de- 
creasing function of halo concentration (§2), the decline of 
Cvir with z implies a relative increase in disc sizes at high 
redshift. This should result in lower than expected surface 
brightnesses at high z, both because of the extended size 
and the corresponding lower efficiency of star formation. To 
this one could add the fact that the supply of cold gas 
for disc formation at high redshifts may be limited (not 
extending all the way to -Rvir) because the smaller inner 
densities will lessen coUisional cooling. These results may 
hinder the association of quiescently star-forming objects 
with Lyman-break galaxies as discussed e.g., by Mo, Mao 



& Whi te (1999 ) (see Kolatt et al. 199^ ; jomcrvillc, Primack 



& Faber 2000 for an alternative physical model for Lyman- 



break galaxies). They further argue for the slow evolution 
of the TuUy-Fisher relation. 



• Disc galaxies at high redshifts are predicted to be more 
extended and of lower surface brightness than expected pre- 
viously. The constant inner radius at fixed mass may be 
a dynamical clue for understanding the formation of halo 
structure. 

• The scatter in the halo mass-velocity relation is signif- 
icantly larger than in the observed TF relation, which sug- 
gests that the luminosity of a disc forming inside a halo of 
a given mass should correlate with the maximum rotation 
velocity. We pointed out a possible simple explanation for 
that. 

• The environmental trends of halo profiles may caution 
against the universality of the TF relation. In addition, these 
trends, together with the observed scatter, may provide in- 
sight into the origin of the Hubble sequence. Below, we ar- 
gue that haloes of low concentration will tend to host blue 
galaxies and haloes of high concentration, red galaxies or 
spheroids. We also point out that extremely low Cvir haloes 
plausibly host LSB galaxies. 



We have proposed an alternative to the toy model orig- 
inally proposed by NFW96. It reproduces the correlations 
between the two parameters of the halo profiles, e.g, Cvir 
and Mvir, as well as the redshift dependence of these cor- 
relations. This model also offers a simple parameterization 
that reproduces the scatter about the median relation ob- 
served in our simulation. The modified toy model is a useful 
tool for semi-analytic modelling of galaxy fo rmation. In par 
ticular, analyses of the type performed by 



Mo et al. 199£ 



which made predictions for disc properties at z ~ 3, based 
on the halo toy model of NFW97, should be reconsidered 
using our modified toy model. A program that implements 
the model for several cosmologies and provides Cvir(Mvir,«) 
with the expected scatter is available from the authors upon 
request. 
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Figure 12. The inner radius Ts as a function of redsliift for dis- 
tinct haloes of fixed mass, Mvir = 0.5-1.0xl0i2;j-iA^g Shown 
are the median (solid line) and intrinsic 68% spread (dashed 
lines). The median value for rs remains approximately constant 
as a function of redshift. 
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Figure 13. The evolution of the distinct halo TuUy-Fisher rela- 
tion, Mvir versus Vmax, for several redshift steps. The evolution 
is weaker than the Mvir versus Vvir relation because Cvir falls 
rapidly with redshift (Fig. pj| ). 

The large intrinsic scatter we find in the correlations 
between the halo profile parameters makes the haloes a two- 
parameter family, as expected, and should be taken into ac- 
count when trying to model the scatter in observable proper- 
ties of galaxies (c.f. Mo, Mao, & White 1998; Navarro 1998; 
Navarro & Steinmetz 1999). 

The spread in exponential disc sizes implied from our 
68% spread in concentration values for an i?vir ~ 200 /i^^kpc 
halo at 2 = (A'/vir ~ 10^^/i"^Mq , Cvir:9.2 21.3) is 
rd : 4.0 2.6 /i~'^kpc (see Eq. This is roughly the same 
spread in disc sizes resulting from a spin parameter variation 



of A : 0.05 0.03, which is approximately th e intrinsic 
spread in A inf e rred from N-body simulations ( Barnes & 



Efstathiou 1987 



Warren et al. 1992 



Bullock et al. 2000a) 




The two quantities, Cvir and A, are thus of comparable 
importance for determining observable properties of galax- 
ies, and define a plane in parameter space for haloes of fixed 
mass. We suggest that the Cvir-A plane can perhaps be linked 
to the observed variations of galaxies and help explain the 
Hubble sequence. For example, we argued above that haloes 
with very low concentrations would tend to lead to discs 
of low surface brightness, and with slowly rising rotation 
curves (i.e., Knax/Kir ^ 1). This argument will only ap- 
ply, however, if A is large enough to prevent excessive gas 
infall to the center of the halo, which would tend to in- 
crease the effective concentration of the system. Similarly, 
haloes of high-Cvir and low-A will likely be unable to pro- 
duce dynamically stable discs (see, e.g.. Mo et al. 1998), 
and instead host spheroids. Other combinations of these pa- 
rameters, may, perhaps, lead naturally to a range of galaxy 
morphologies. 

This kind of mapping is further motivated by the under- 
standing that high-Cvir haloes collapse earlier than low-Cvir 
haloes (as predicted by the toy model and explicitly demon- 
strated by NFW97 using simulated haloes). A natural as- 
sociation is then that high-Cvir haloes host old, red galax- 
ies, and lower-Cvir haloes host young, blue galaxies. Further- 
more, the environmental trend, that haloes in low-density 
environments tend to be less concentrated than haloes of 
the same mass in high-density environments, fits nicely into 
this picture. Indeed, LSB galaxies are observed to be more 
isolated than galaxies of hi gher surface brightness (Bothun 



Mo et al. 1994 ), and spheroids tend to inhabit 



high- density environments (Dressier 198C ; Postman & Geller 



We have made in this paper only crude preliminary at- 
tempts to study the implications of our results concerning 
halo profiles. In Bullock et al. (2000) we have performed 
a first step towards understanding how Cvir evolution af- 
fects the nature of the galaxy luminosity-velocity relation 
at high redshift. It would be desirable to go on and incor- 
porate our measured halo properties, including scatter, into 
semi-analytic modelling of gas processes and star formation 
in order to make detailed predictions for observable galaxy 
properties and their evolution. 
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APPENDIX A: THE NFW97 MODEL 

For completeness, we briefly review the NFW prescription 
for determining Cvir{M, a). The goal of the NFW procedure 
is to provide the density profile of a halo of mass Mvir at 
epoch ao assuming that the profile is of the NFW form (|l|). 
The collapse redshift is now determined using the Press- 
Schechter approximation, which, given Mvir at epoch ao, 
can be used to approximate the probability distribution for 
the epoch a when a halo trajectory was first more massive 
than some fraction / of Mvir (Lacey & Cole 93) 



P(> /Mvir, a|Mvir, ao) = erfc 



'?crit(a) — Scritjl) 
V2(a2(/Mvir)-a2(Mvir)) 



Uc, by solving Eq. 
setting P(> /M, 



The quantities ao and 5crit are defined in § y. One determines 
Al for the most probable value of a by 
tc|M, ao) = 0.5. One now assumes that 
the central density of the halo is proportional to the density 
of the universe at Uc, which implies: 



the virial radius of large haloes, therefore certain ambigui- 
ties arise. How close must two density maxima be in order 
for them to represent a single object? How does one differ- 
entiate substructure from a collision in progress? How does 
one assign mass to haloes and subhaloes appropriately? 

We have found a solution to these problems by assign- 
ing to each halo two length scales — an inner radius, r^, and 
and an outer radius _Rvir . We do so by modelling the density 
profile of each halo using Eq. ^. The virial radius JJvir deter- 
mines each halo mass and radial extent, and rs determines 
when two density maxima/haloes should be combined into 
one. The details are described below. 
.(Al) Because the modelling process requires fitting a density 
profile to each halo, we attempt to find only haloes with 
more than A'^^"" particles within their virial radii. If rUp is 
the mass of each particle, this means the minimal virial mass 



Pb = kpu{ao) 



ao 



(A2) 



where fc is a numerical constant. Now, given Mvir and ao, 
and A.2 determine ps and thus completely specify 



Eqs. Al 



the density profile. 

This procedure has two free parameters, / and k, which 
may be adjusted in order to match the slope and normal- 
ization respectively of Cvir(M) at a = 1. NFW97 show that 
this two parameter model reproduces the a = 1 relations 
of simulated haloes in various cosmologies, including power 
law and open models. For the ACDM model we discuss, 
their favorite parameters are / — 0.01 and k = 3.4 x 10'', 
and these provide a reasonable reproduction of the median 
Cvir(Af ) relationship at a = 1 in our simulations. 

Although useful in its ability to provide the correct 
2 = relation, this model fails to reproduce the observed z 
evolution of halo concentrations (§ ^). If the value of / is ad- 
justed to the unphysically small value of / = 10^^^ then the 
redshift behavior is becomes similar to that reported in § ^ 
However, this value of / pushes the relevant region of the 
power spectrum to extremely small masses, where the effec- 
tive slope is extremely fiat, as0.05, so the implied Cvir(M 
dependence becomes shallower than what is observed. In 
we present a revised toy model which, using the same num- 
ber of free parameters, reproduces the full observed behavior 

of Cvir(M, Z). 



APPENDIX B: THE HALO FINDER 

Most commonly used halo finders, which work either by the 
location of overdensities in a spatial window of fixed shape 
(usually spherical) or by friends-of-friends algorithms, do not 
account for haloes within haloes. Since our projects specif- 
ically address the question of substructure, we have been 
obliged to devise a halo finder and classification algorithm 
suited for this purpose. 

If one were only interested in distinct, virialized objects, 
haloes would be easily identified — there is little confusion as 
to where one halo ends and another begins because the phys- 
ical extent of an object is determined by the virial overden- 
sity criterion. However, we are interested in objects within 



of haloes identified is M^^'i" 



Equivalently, 



using Eq. |2|, we have a minimum virial radius -R™r". The 
value of N^^'^ is the first free parameter of this algorithm. 
We use A^;^"" = 50. 

O ur density maxima findin g routine is based on the 
BDM (|Klypin fc Holtzman 1997|) algorithm. We outline our 



procedure below, including our precise methodology which 
is generalizable for any simulation parameters and our de- 
tailed procedure for defining and classifying haloes. 

(i) We constru ct density field values by a Cloud-in-Cell 
(CIC) process ( Hockney fc Eastwood 1981 ) on the largest 
grid of the simulation AL, and rank the particles according 
to their local density as determined on this grid. 

We then search for the possible halo centres, using two 
sets of smoothing spheres; one, with a small radius, rspi, in 
order to locate of tight, small clumps; and the other, with 
a larger radius, rsp2, in order to locate the centres of haloes 
with diffuse cores. The smaller radius is rspi = a /res, where 
/res is the highest force resolution in the simulation and a is 
a free parameter of order unity. We use a = 2. The second 
set of spheres have rsp2 ~ ^"ir"- 

For each set of spheres, we take from the ranked list the 
particle with the highest local density and place a sphere 
about its location. A second sphere is placed about the next 
particle in the list not contained in the first sphere. The 
process is continued until all of particles are contained within 
at least one sphere. Because we are only interested in centres 
of haloes more massive than M™^", we discard each sphere 
with fewer than a set number of particles. The minimum 
number of particles required for a kept sphere is determined 
separately for each radius. 

For the r-spi spheres, we use the following conservative 
halo density profile: 



p{r) 



C/ ''spl 



r < rspi 
r > rspi. 



(Bl) 



(where C is determined my fixing the minimum halo mass 
to be M^^.'^), in order to estimate the minimum core number 
of particles within rspi: 



A"spi 



JVp 



l + 6[(i?™A=p)^/^-l]' 



(B2) 



For the z — timestep of the 60/i ^Mpc simulation we an- 
alyze, Aspi = 3.9 3. Spheres of size rspi with fewer than 
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Nspi particles are discarded. Similarly, all of the rsp2 spheres 
containing fewer than A'gp2 = A^™'" particles are discarded. 

The final list of candidate halo centres is made up of all 
of the (small) rspi spheres, together with each of the rsp2 
spheres that does not contain an rspi sphere. 

(ii) For each sphere of radius r^p = rspi or rsp2, whichever 
applies, we use the particle distribution to find the centre of 
mass and iterate until convergence. We repeat the proce- 
dure using a smaller radius, r — Vi, where = We 
continue this method until ri = r^, where tl is defined by 
the criterion r_L > /res > fL+i, or until reduction leads to 
an empty sphere. 

(iii) We unify the spheres whose centres are within 

of each other. The unification is performed by making a 
density weighted guess for a common centre of mass, and 
then iterating to find a centre of mass for the unified object 
by counting particles. The size of sphere used to determine 
the centre of mass is the smallest radius that will allow the 
new sphere to entirely contain both candidate halo spheres. 

(iv) For each candidate halo centre we step out in radial 
shells of 1 h~^kpc, counting enclosed particles, in order to 
find the outer radius of the halo: Rh = min(_Rvir, Rt)- The 
radius R^ir is the virial radius, and Rt is a "truncation" 
radius, defined as the radius (< -Rvir) in which a rise in 
(spherical) density is detected (dlogp/dlogr > 0). This is 
our method for estimating when a different halo starts to 
overlap with the current halo and is important for haloes in 
crowded regions. We estimate the significance of a measured 
upturn using the Poisson noise associated with the number 
of particles in the radial bins considered. Only if the signal 
to noise of the upturn is larger than an^ do we define a 
truncation radius. The value of cthj is a free parameter. We 
use (Tflt = 1.5. 1^ 

(v) Among the halo candidates for which we have found 
an Rvii, we discard those with A'^vir < A^™'", where A'vir is 
the number of particles within -Rvir. Among the halo candi- 
dates for which we have found a rise in spherical density, we 
discard those which contain less than A^^'" particles, where 

J^rn^„ ^ ^min ^ ^min^ otherwise 



7» rmzn ■p.rrmn 

J^Rt - J^p 



Rt 

Dmin 



(B3) 



The above constraint follows from an extrapolation of the 
minimum mass halo using an isothermal profile p{r) oc 

(vi) We model the density profile of each halo using 
the NFW form (Eq |l|) and determine the best fit 
and pa values, which determine -Rvir and -Mvir. The fit- 
ting procedure uses logarithmically spaced radial bins from 
max(2/res, 0.02 x mm{Rvii, Rt)) out to Rh- If any bins are 
empty we decrease the number of bins by one until all bins 
are full. If the number of bins is reduced below 3 we discard 
the halo as a local perturbation. 

The fit takes into account the Poisson error in each bin 
due to the finite number of particles, and we obtain errors 
on the fit parameters {ar, and (jp^) using the covariance 



M The choice was motivated by several tests using mock cata- 
logues of haloes in clusters designed to determine how varying 
fJijj affects our ability to fit the density profiles of subhaloes. Al- 
though our results were not strongly dependent on this choice, we 
did obtain the best fits using an^ =1.5. 



matrix in the fit routine. The errors on the fit parameters 
can be translated easily into errors for -Rvir (f-Rvir) ^^'^ the 
estimated NFW mass of each halo, A/vir (ctm). 

(vii) We unify haloes which overlap in R^. Our criterion is 
met if two (or more) halo centres have R^ radii which overlap 
with each other while at the same time having velocities 
which allow them to be bound to the common system. If 
such a case occurs, then along with the individual halo NFW 
fits, we fit another NFW profile about the common centre 
of mass of the two combined haloes and decide whether the 
candidate- united-haloes are bound/unbound to the common 
NFW fit using the radial escape velocity determined using 
the common NFW profile (see below). If both haloes are 
bound we combine the two haloes into one, and keep the 
common fit for the characteristic parameters. If at least one 
is not bound, we do not combine the haloes. 

An exception to this unifying criterion occurs if the fit 
errors on -Rs are large {ur^/Rs > 1), we replace -Rs — > 
min(Rs, Rt). In addition if the -Rs of a halo obeys -Rs > -Rvir" 
then we relax the strict combining of overlapping haloes. 
This case, which we refer to as the "cD" halo case, is dis- 



cussed below (see (ix).) 

(viii) For each halo, we remove all unbound particles be- 
fore we obtain the final fits. We loop over all particles within 
the halo and declare a particle at a distance r from the centre 
of a halo to be unbound if its velocity relative to the centre 
of mass velocity of the halo obeys v > ^2 |<I>nfw(7')| , where 

the radial potential for NFW density profile is given byPj 



$NFw(7-) = -47rGps-Rs 



log(l 



(B4) 



After removal, we construct a new density profile and find 
new NFW fit parameters. The procedure is repeated un- 
til the number of unbound particles becomes < 1% of the 
bound particles or until -M";™ < -M^'i", in which case the 
halo is discarded. 

An exception to this removal scheme occurs if two haloes 
lie within the virial radius of each other and the ratio of their 
meisses are at least 0.75. We define haloes in this situation 
to be a "partner" pair. For each halo in this situation, we 
take not only its potential into account, but also that of its 
partner. 

(ix) An interesting case of subhalo structure, which would 
otherwise be excluded from our finding algorithms, is that of 
one or more density peaks close to the centre of a large halo. 
We shall refer to these inner density peaks as cD cases. If a 
halo, after its unbound particles have been removed, obeys 
the following criteria, it is a candidate for containing cD 
haloes: a) the NFW fit has a standard GoF < 0.1, b) the 
halo is a host of at least one subhalo, and c) the halo is 
"large", with rs > JJ^j.". 

We identify the potential centres of cD haloes by searching 
the bound particle distribution within rs of each candidate 
cD halo using a CIC process on a fine grid (rspi = a/rcs). We 
discard all candidate density peaks with local densities less 
than the extrapolated minimum density (above the back- 



•58 Note that this potential is not necessarily the physical gravita- 
tional potential at the halo location. For a subhalo, for example, 
the host background potential is not included. 
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ground density) within the core region of our smallest halo 
(see item |(i)| .) 

For each density maxima located farther than rspi from 
the centre of the candidate host halo, we find Rt and fit a 
NFW profile with iterative unbound particle removal. These 
are our cD haloes. cD haloes are discarded if their extrapo- 
lated virial mass is lower than our minimal mass halo. 

Because we have a strict mass limit Af™^" — TV™'" x nip, 
we expect our halo finder to be somewhat incomplete just 
above M™^". We have also checked our completeness in two 
ways. First, we used a separate BDM halo finder that does 
not attempt to fit profiles and does not demand the unifi- 
cation of haloes within a specified radius. It does, however, 
unify haloes that have equal velocities withi n 15% as long as 
they have centres within ~ 150 ^~^kpc (see Klypin & Holtz- 
man 19p7|). This procedure allows a complete identification 
of DM haloes down to much lower particle numbers than 
our own. In order to check our results we have assigned to 
each halo in the catalog produced by the other finder a typ- 
ical Vs given its mass, and checked the returned halo list for 
consistency against our catalogue from the same simulation 
box. For A''particics ~ 150, we estimate ~ 80% completeness 
and for A^particics ~ 500 we obtain a95% completeness. A 
second, and almost identical completeness determination is 
obtained by car efully analyzing th e roll-over in our observed 
mass function ( ^igad et al. 2000 ). We attribute our incom- 
pleteness for small masses to our fitting procedure, and er- 
rors associated with this process. 



APPENDIX C: REDUCING TF SCATTER 

As discussed in §6, we find for distinct haloes a la scatter 
at fixed Mvir of AVmax/Knax — 0.12, which is between 1.3 
and 4 times the range of the reported intrinsic TF scatter. If 
these ACDM haloes are to host galaxies like those observed, 
this excessive scatter must be reduced. The translation of the 
halo virial mass into a disc luminosity, and of the original 
halo Knax into a final observed disc velocity, should some- 
how decrease the scatter. Following is a qualitative analysis 
of how this can come about in a natural way. The idea is 
that for a fixed halo mass and spin, a higher Vmax should 
induce further gas contraction into smaller radii, and there- 
fore higher gas density, star-formation rate and luminosity. 
This can be shown in a little more detail, as follows. 

The size of the exponential disc, r^, that forms by a dis- 
sipative contraction of gas inside a given dark-matter halo 
can be estimated under the adiabatic baryonic-infall approx- 
imation. We showed in §2 (see Equation ^ that ra is a de- 
creasing function of Cvir for haloes of a fixed virial mass and 
spin, as long as the disc mass is a constant fraction of Mvir. 
For a typical case of a Vvir ~ 200 km/s halo, we demon- 
strated that, in the range encompassing 68% of Cvir for such 
haloes (cvir:9.2 — > 21.3), the corresponding spread in disc 
sizes is Rd : 4.0 -+ 2.6 h~^kpc. Across this range, an effec- 
tive power-law approximation would therefore be rd cx c~;^'^. 
If (a) the gas density in the disc scales like p oc r^^, (b) the 
star formation rate obeys a typical Schmidt law, p oc p^'^, 
and (c) the luminosity scales like L oc pr^, then the lumi- 
nosity at a given mass depends on Cvir as L oc r^^ oc c^\^. 

is also a monotonic function 



Since Vmax 

of Cvir (Eq. 



'for a fixed Vvir, 
1; the effective power-law approximation across 



the 68% range is Vmax oc Cvir ), we have obtained a positive 
correlation between L/Mvir and the deviation of Knax from 
the median Vmax at a given Mvir. Ignoring, for the moment, 
any difference between the Vmax of the original halo and 
that of the disc, the obtained correlation would correspond 
to a reduced scatter in L by a factor larger than 2, as re- 
quired. The effect of the spread in spins on the TF scatter is 
expected to be reduced for similar reasons, namely, because 
the luminosity and the maximum velocity are both expected 
to correlate with the spin in the same sense. 

More detailed modelling, which takes into account how 
Vinax changes as the baryons fall in, will be needed to test 
this hypothesis in detail. 
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